################################################################## 
## 全局变量

export work_dir=/public/home/xxf2019/20230331_colrNew/public_data/liver
export dna_preprocess_scripts=~/dna_preprocess_scripts
export config_path=${work_dir}/config
export config_file=${config_path}/config.sh

source ${config_path}/config.sh

################################################################## 
## contamination asessment
rm -rf ${qc_path}/Summary_contamination.table
for file in ` ls ${vcf_path} | grep con `
do
echo $file
cat ${vcf_path}/${file} | grep -v sample | awk '{OFS="\t"}{print $1,$2}' >> ${qc_path}/Summary_contamination.table
done
sed -i '1i Sample\tContamination' ${qc_path}/Summary_contamination.table

################################################################## 
## 统计所有样本的mosdepth计算的平均深度
echo -e "Tumor\tNormal\tTumorDepth\tNormalDepth" > ${qc_path}/Depth_summary.txt
for line in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo ${line}
Tumor=`echo ${line} | awk -F',' '{print $1}' `
Normal=`echo ${line} | awk -F',' '{print $2}' `
TumorDepth=`cat ${qc_path}/bam/mosdepth_QC/${Tumor}/${Tumor}.mosdepth.summary.txt | grep total_region | awk -F'\t' '{print $4}'`
NormalDepth=`cat ${qc_path}/bam/mosdepth_QC/${Normal}/${Normal}.mosdepth.summary.txt | grep total_region | awk -F'\t' '{print $4}'`
echo -e "${Tumor}\t${Normal}\t${TumorDepth}\t${NormalDepth}" >> ${qc_path}/Depth_summary.txt
done

################################################################## 
## 计算覆盖超过10x的区域
## 为计算突变负荷
## 最后的突变只看了1-22号染色体上的，因此只计算1-22号染色体的区域覆盖10x的区域长度
echo -e "Tumor\tcoverage" > ${qc_path}/Burden.coverage10x.Autosomal.txt

for sample in `cat ${config_path}/tumor_normal.list  | grep -v Tumor | awk -F',' '{print $1}'` 
do
echo $sample
region=`zcat ${qc_path}/bam/mosdepth_QC/${sample}/${sample}.thresholds.bed.gz | awk -F'\t' '{if($1!="X" && $1!="Y")print}' | awk -F'\t' '{sum+=$6}END{print sum}'`
echo -e ${sample}"\t"${region} >> ${qc_path}/Burden.coverage10x.Autosomal.txt
done

################################################################## 
## CDS区域的覆盖深度
## 只纳入常染色体
export region_type=cds
export capture_bed=${ref}/PCAWG_Elements/web_hg19/gc19_pc.cds.use.nochr.bed

for sample in `cat ${config_path}/tumor_normal.list  | grep -v Tumor | awk -F',' '{print $1}'` 
do
echo ${sample}
sh ${dna_preprocess_scripts}/mosdepth_run_region.sh ${sample} ${config_path} ${capture_bed} ${region_type}
done

## 合并cds深度到一个文件
echo -e "Tumor\tcoverage" > ${qc_path}/Burden.coverage10x.Autosomal.cds.txt
for sample in `cat ${config_path}/tumor_normal.list  | grep -v Tumor | awk -F',' '{print $1}'` 
do
echo $sample
region=`zcat ${qc_path}/mosdepth_cds/${sample}.thresholds.bed.gz | awk -F'\t' '{if($1!="X" && $1!="Y")print}' | awk -F'\t' '{sum+=$6}END{print sum}'`
echo -e ${sample}"\t"${region} >> ${qc_path}/Burden.coverage10x.Autosomal.cds.txt
done


################################################################## 
## GGA
## 需跑完vcf质控
conda activate test_new

for Normal in ` cat ${config_path}/tumor_normal.list |  grep -v Normal | awk -F, '{print $2}'  | sort -u ` 
do
echo $Normal
sh ${dna_preprocess_scripts}/GGA_1_CombineVcf.sh ${Normal} ${config_path}
done

## Recal
for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `
echo " sh ${dna_preprocess_scripts}/GGA_2_Recall.sh ${Tumor} ${Normal} ${config_path} " | \
qsub -N ${Tumor}_"GGA" -l nodes=1:ppn=10,mem=10gb,walltime=240:00:00 -q batch -d ${Qsub_log_path}
done

## GGA 只保留原来的突变，去除多等位原来未通过的位点
echo Sample","Raw_Num","GGA_Abstract_Num > ${qc_path}/GGA_Filter.list
for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `

echo "sh ${dna_preprocess_scripts}/GGA_3_Filter.sh ${Tumor} ${Normal} ${config_path}" | \
qsub -N ${Tumor}_"GGA_Filter" -l nodes=1:ppn=2,mem=10gb,walltime=240:00:00 -q student -d ${Qsub_log_path}
done

## Functator注释 GGA
cat ${config_path}/tumor_normal.list | grep -v Normal | awk -F, '{print $1,$2}' | xargs -P 10 -I sample sh -c '
sh ${dna_preprocess_scripts}/GGA_4_Functator.sh sample ${config_path}
'

## 合并同一人的vcf，用于进化树分析
muti_sample=`cat ${config_path}/tumor_normal.list | awk -F, '{print $2}' | sort | uniq -c  | awk '{if($1>1){print $2}}' | tr '\n' '|' | sed 's/|$//'`
echo "Normal,Raw_mut_num,Tree_mut_num" > ${qc_path}/Treeomics_Mut_Num.list
cat ${config_path}/tumor_normal.list | grep -E ${muti_sample} | grep -v Normal | awk -F, '{print $2}' | sort -u | xargs -P 5 -I Normal sh -c '
echo Normal
sh ${dna_preprocess_scripts}/Treeomics_CombineVcf_v4.sh Normal ${config_path}
'

################################################################## 
## 统计纯度和CNV
rm -rf ${qc_path}/Purity_facets.csv
rm -rf ${facets_path}/Facets_all_seg.seg

for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}'  `
Normal=` echo $samples | awk -F, '{print$2}' `
sh ${dna_preprocess_scripts}/Facets_PuritySummary.sh ${Tumor} ${Normal}
done

sed -i '1i Sample,Purity,Ploidy' ${qc_path}/Purity_facets.csv
sed -i '1i ID\tchrom\tloc.start\tloc.end\tnum.mark\tseg.mean' ${facets_path}/Facets_all.seg


################################################################## 
## MSI
cat ${config_path}/tumor_normal.list | grep -v Normal | awk -F, '{print $1,$2}' | xargs -P 5 -I sample sh -c '
sh ${dna_preprocess_scripts}/Msisensor.sh sample 
'

## 整理总结输出
rm -rf ${qc_path}/All_Msi.tsv
for samples in ` cat ${config_path}/tumor_normal.list | grep -v Normal | tr '\t' ',' ` 
do
echo $samples
Tumor=` echo $samples | awk -F, '{print$1}' `
Normal=` echo $samples | awk -F, '{print$2}' `
cat ${MSI_path}/${Tumor}_${Normal} | sed '1d' | awk '{print Normal"\t"Tumor"\t"$0}' Tumor=${Tumor} Normal=${Normal} >> ${qc_path}/All_Msi.tsv
done
sed -i '1i Normal\tTumor\tTotal_Number_of_Sites\tNumber_of_Somatic_Sites\tRate' ${qc_path}/All_Msi.tsv
